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ABSTRACT 

We model strong forward shocks in young supernova remnants with efficient particle acceleration 
where a nonresonant instability driven by the cosmic ray current amplifies magnetic turbulence in the 
shock precursor. Particle injection, magnetic field amplification (MFA) and the nonlinear feedback of 
particles and fields on the bulk flow are derived consistently. The shock structure depends critically on 
the efficiency of turbulence cascading. If cascading is suppressed, MFA is strong, the shock precursor 
is stratified, and the turbulence spectrum contains several discrete peaks. These peaks, as well as the 
amount of MFA, should influence synchrotron X-rays, allowing observational tests of cascading and 
other assumptions intrinsic to the nonlinear model of nonresonant wave growth. 
Subject headings: acceleration of particles - shock waves - cosmic rays - supernova remnants — 
magnetic field - instabilities 



1. INTRODUCTION 

Diffusive shock acceleration (DSA) in supernova rem- 
nant (SNR) shocks is the most likely mechanism for the 
origin of the bulk of galactic cosmic rays (CRs) up to 
at least 10^^ eV. To produce the observed density of 
galactic CRs, the acceleration mec hanism must be effi- 
cient and, therefore, nonlinear (e.g..lBlandford fc EichleJ 
[l987t iMalkov fc Drurvl [20011 : iHillasi 12005') . X-rav imag- 
ing of young SNRs reveals thin synchrotron filaments 
that are best explained as originating from shock ac- 
celerated TeV electrons radiating in strongly amplified 
circumstellar magnetic fields fsee [V ink & Lamins 200l 
iBamba et aI1l2003l:IUchivama et all2007l: , Vink 2008, and 
references therein) . These fields are almost certainly pro- 
duced as part of the DSA proce s s (e.g . . iBell &: Lucek 
Ml iBelll [200l fAmato fc Blasil [200a iPelletier et al " 
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iBell ()2004D suggested that magnetic fluctuations in a 
shock precursor may be amplifled by a fast, nonresonant 
instability induced by the strong diffusive electric cur- 
rent of the shock-accelerated CRs. The transverse MHD 
waves produced by this mechanism have wavelengths 
much shorter than the wave generating CR gyroradii. In 
efficient DSA, the CR current driving the instability must 
be strongly coupled to the turbulent magnetic fields and 
supersonic MHD fiow. To address this issue, we simu- 
lated the amplification of Bell's modes as^part of a nonlin- 
ear Monte Carlo (MC) model of DSA (jVladimirov et all 
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|2006[ ). We self-consistently model four basic strongly 
coupled components of the system: the bulk plasma fiow, 
the full particle spectrum, the self-generated MHD tur- 
bulence including cascading, and thermal particle injec- 
tion. The efficiently produced CRs modify the upstream 
flow speed, and the CR current determines the growth 
of stochastic magnetic fields. The fields, in turn, set the 
particle diffusive transport properties and, subsequently, 
the injection and acceleration efficiency of the particles, 
closing the system. 

2. THE MODEL 

Our Monte Carlo model of nonlinear particle accel- 
eration with strong, r e sonan t, MFA is presented in 
(jVladimirov et al.ll2006l [20081 ). Here, we replace reso- 
nant growth with the fast nonresonant instability of Bell 
and include the cascading of that turbulence to shorter 
wavelengths. The momentum and space dependent par- 
ticle mean free paths are calculated consistently with the 
turbulence. 

Consider the nonlinear precursor of a strong, plane- 
parallel, steady-state shock. In the reference frame co- 
moving with the subshock (at a; = 0), where the plasma 
flows in the positive a:-direction, the flow speed u{x) 
has the value uq far upstream. The speed, u(x), drops 
in the precursor until it obtains the downstream speed 
U2 — uo/-Rtot, where -Rtot is the overall shock compres- 
sion ratio. The nonresonant current instability is as- 
sumed to be the source of the wave spectrum energy 
density W(x, k) (here k is the wavenumber, and Wdk 
is the amount of energy in the interval dk per unit spa- 
tial volume) according to the following equation: 
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The quantity Tl{x, k) describes cascading, i.e., the trans- 
fer of turbulence energy from long to short wavelengths, 
L(x, k) represents the dissipation of turbulence, and 
7nr(a;, k) is the quasi-linear rate of wave ene rgy amplific a- 
tion by the instability, and is given by (see iBelll (|2004J l). 



7nr = 2vAk\ 7^-1, for 1/rgi < k < kc- (2) 



Here VAix) — Bq/ y^ATTp{x) is the Alfven speed, c is the 
speed of light, Bq is the far upstream magnetic field di- 
rected towards the shock normal, p{x) is the thermal 
plasma mass density, rgi{x) is the gyroradius of the least 
energetic current generating CR, the critical wavenum- 
ber kc{x) — Airjdix) /{cBq), and the local diffusive elec- 
tric current of CRs responsible for the instability, jd{x), 
is determined by the MC simulation. 

Bell's derivation of ^ assumes that the gyroradii of 
the streaming CRs are much greater than the wave- 
lengths of the generated waves, which is expressed here 
by the condition 1/rgi < k. The growth rate, 7nr, has a 
maximum at fc = kc/2 and vanishes for k > kc- 

Nonlinear interactions between turbulent harmonics 
may lead to a re-distribution of the energy of turbulent 
fluctuations in fc-space (i.e., spectral energy transfer). 
The term dli/dk in Equation ([1]), with 11 > 0, describes 
a transfer of energy from large to small turbulent fluc- 
tuations (i.e., cascading). Such a description, i.e., the 
Kolmogorov cascade model, has been successful in ex- 
plaining the spectra of locally isotropic, incompressible 
turbulence in non-conducti ng fluids, observ ed in experi- 
ments and simulations (e.g. . lBiskamp|l2003l ) . However, it 
is uncertain how spectral transfer operates in a collision- 
less shock precursor with a CR current strong enough to 
modify the MHD modes, and in the presence of strong 
magnetic turbulence. In weak MH D turbulence, cascad- 
ing w as shown to be anisotropic ()Goldreich fc Sridhail 
Il995l l: harmonics with wavenumbers transverse to the 
mean magnetic field experience a Kolmogorov-like cas- 
cade, while the cascade in wavenumbers parallel to the 
mean field is suppressed. [Diamond fc Mal kov (2003) ar- 
gue for an inverse cascade (i.e., from short to long wave- 
lengths) of turbulent harmonics in the context of mag- 
netic field amplification in DSA. Spectral energy transfer 
in the form of anisotropic d iffusion of energy in fc-space 
was proposed by Matthaeu s et al.l (pOOS) . MHD model- 
ing ofBeirsjnstabilit;Y_w2tlifi^ performed 
bv iBeU (|2005| ) and IZirakashvili et all jlOOl), show that 
the magnetic structures consist of expanding walls of 
strong magnetic field spiraling around cavities in density 
and magnetic field. Moreover, Zirakashvili et al. (2008) 
revealed that the turbulent energy is non-linearly trans- 
ferred to both longer and shorter scales (see their Figure 
4). It should be noted however, that the MHD mod- 
els ignored both the interaction of CR particles with 
the instability-generated turbulence, and the feedback 
of the turbulence on the CR current that are the main 
subjects of our study. The deflection of CRs by mag- 
netic fluctuations may lead to important nonlinea r ef- 
fects, such as a saturation of the instability (e.g., iBelll 
l2005t iPelletier et al.l l2006f l . A complete numerical de- 
scription of turbulence, that includes the evolution of 
the fields and the motion of the particles, with a realis- 
tic dynamic range, is currently unfeasible. Because the 



primary goal of our present research is to study the non- 
linear aspects of DSA with Bell's instability including the 
interaction of CR particles with the instability-generated 
turbulence, approximations must be made for the poorly 
known details of the turbulent spectral energy transfer. 
In this work, we investigate the consequences of different 
spectral energy transfer regimes by presenting two limit- 
ing cases. In one (shown with solid curves in all plots), 
cascading is fully suppressed, i.e., H = 0. In the other 
(shown with dotted curves in all plots), the cascading 
from large to small scales is efficient and has the differ- 
ential form (see, e.g.. IVermall2004D . H = W^3/2^5/2^-i/2^ 
corresponding to Kolmogorov's model^. 

The dissipation term, L, is assumed to be zero with 
no cascading and to have the form of viscous dissipa- 
tion, L = v^fc^fcT^VF, in the model with cascading (e.g., 
IVainshtein et al.ril993f ). The wavenumber, kd, is iden- 
tified with the inverse of a thermal proton gyroradius: 
kd(x) = eBQl[c^ mpkBT{x)] , where nip is the proton 
mass, fcs is Boltzmann's constant and T(x) is the local 
gas temperature d etermined from th e Jieating induced by 
L, as described in I Vladimirov et all (|2008 h. We will con- 
sider elsewhere models of spectral transfer different from 
Kolmogorov cascading. 

To highlight the effects of the Bell nonresonant insta- 
bility and cascading in strong SNR shocks of SNRs, we 
here ignore reson ant instabilities and th e compression 
of turbulence (see IVladimirov et al.ll2006l ). even though 
thes e effects may play an important role in weaker shocks 
(see lPelletier et al1 l2006V We have performed more com- 
plete calculations including resonant and nonresonant in- 
stabilities with compression and cascading (to be demon- 
strated elsewhere) and confirmed that the qualitatively 
new findings presented here stand out even when reso- 
nant instabilities and compression are included. 

Equation ([T]) is integrated with respect to x using the 
far upstream (x = — oo) boundary condition W oc k~^. 
We assume that these far upstream seed fiuctuations are 
linear waves that are not subject to cascading or dissi- 
pation, and that the transition to the turbulent regime 
takes place at a position xt where the amplified wave 
spectrum reaches the value kW{xt,k) = Bq/Att at some 
fc. At this point, in the model with cascading, H and L 
are set from zero to the values defined above. 

The mean free path of a particle with momentum 
p, A(x,p), is determined by the MHD turbulence. We 
calculate A by combining well known theories in differ- 
ent parameter ranges. For the highest energy particles 
with gyroradii, r^, that exceed the turbulence correla- 
tion length IcoY, we assume A = r^//cor oc p^, as de- 
scribed in iToptvginI (|1985D (see also lJokipiil (|197lD '). Par- 
ticles with rg that resonate with the field fluctuations are 
assumed to have A — TgjT (e.g., iBlandford fc Eichleil 
Il987f ). where T is the normalized product kW evaluated 
at fc resonant with the particle momentum p. The low- 
est energy particles, with Vg <C ^con are assumed to be 
trapped by turbulent vortices, and their diffusion is de- 
termined by the dynamics of these vortices; these parti- 
cles are assumed to have A = Zcor, independent of p (e.g., 
iPalmeil ()1982f )). Our model for A smoothly interpolates 
between these regimes by introducing a wavenumber fc* 

^ Note that we ignore a c onstant on the order of unity in this 
expression fe.e.. IVermall20(3 ). 
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that, for a given particle momentum p, separates the tur- 
bulence spectrum W into the small-scale (fc > /c,) and 
the large-scale {k < k^,) parts. In addition, we assume 
that thermal particles not yet injected into the accelera- 
tion process have the Bohm mea n free path A = Tq . Full 
details of the model are given in IVladimirovl (l2009h . 

The diffusive CR current jd(x) is calculated as the cor- 
responding moment of the accelerated proton distribu- 
tion function f{x, p) derived by the MC simulation. 

At every upstream location x, the minimal CR momen- 
tum pi is defined from the particle distribution function, 
/(a;,p), as the momentum below which the CRs con- 
tribute 1% of the total CR pressure. Then pi is used to 
calculate the gyroradius rgi that determines the applica- 
bility range of ([2]). We assume 7,11. = outside of that 
range of k. 

3. RESULTS 

For both examples shown here, we model a shock of 
speed uq = 10* km s~^ propagating in a uniform mag- 
netic field i?o = 3 fiG in a plasma with proton density 
no — 0.3 cm"'^ and temperature Tq — 10"* K. We assumed 
that the acceleration process is in steady-state and size- 
limited with a free escape boundary at Xfcb = — lO^fgo, 
where Vgo = muQc/eBo w 3.5 x 10^" cm. Note that Pmax 
scales with Xfeb- In an actual shock, the finite size, 
and/or finite age, of the shock will determine Pmax- 
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Fig. 1. — Flow speed u, effective magnetic field -Bcft , diffusive 
CR current j^, and temperature T vs. position x. SoHd curves are 
without cascading and dotted curves include strong Kolmogorov 
cascading. The viscous subshock is at x = and the position of 
the free escape boundary (FEB) is indicated. 

In Fig. [T] we compare the flow speed, u{x), the ampli- 
fied effective field, Bcs{x), the CR current, jd(x), and 



the thermal plasma temperature, T{x), for the two mod- 
els. With no cascading (solid curves), the overall com- 
pression ratio i?tot = uo/u2 ~ 15, and the downstream 
field, Bcs{x > 0) « 1000 /xG. These quantities are sub- 
stantially reduced with cascading (dotted curves), i.e., 
^tot ~ 11, and Bes{x > 0) ~ 100 //G. There are also dif- 
ferences in jd{x), where jd is smooth with cascading but 
shows an uneven structure without cascading. The pa- 
rameter that determines the scale of the generated turbu- 
lence, fee, has the same x-dependence as jd- This uneven 
precursor structure has a striking effect on the magnetic 
turbulence (see Fig. [3]). Another prominent difference is 
the significantly increased precursor temperature, T(a:), 
with cascad ing due to the dissipat ion of turbulence at 
large k (see IVladimirov et all (|2008( )). 
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Fig. 2. — Proton distribution function, /, times p*, and accel- 
erated particle mean free path, A, both calculated downstream of 
the subshock. The turnover momenta, Pmax, for the two cascading 
models are indicated. 

In Fig. [2] we plot the downstream (x = 5rgo) parti- 
cle distribution p*/(x,p) and mean free path X(x,p)/rgQ. 
The maximum momentun, Pmax, is clearly greater with- 
out cascading (~ 10^ vs. '--^ 2 x 10"* nipc) but the con- 
cave shape for p"*/ is present in both cases, indicat- 
ing that both are effici ent particle accelerators (e.g., 
iBerezhko fc Ellis"oiilll999f l. As indicated by the thermal 
peaks, the shocked temperature is considerably higher 
in the model with cascadi ng, an e ffect similar to what 
was observed in Iviadimirov et all (|2008[ ). With cascad- 
ing (dotted curves), A is a smooth function of p showing 
the three regions discussed above: A ~ constant for low 
p, X oc p for intermediate values of p, and A oc for the 
highest momenta. Without cascading, A is always less 
than with cascading and shows plateau regions where 
particles are trapped by vortices of different scales. 
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Fig. 3. — Turbulence spectra, W , times k, at different locations 
relative to the subshock. The lightest- weight curves are calculated 
at Xfob = — fO^ TgQ, the heaviest curves are calculated in the down- 
stream region, and the medium-weight curves are calculated in the 
precursor at x = — 1.2x fO* VgQ. The far upstream seed turbulence 
(at all k) is at the level indicated by the horizontal lines. 

4. DfSCUSSfON 

Our most intriguing result is shown in Fig. [3] where 
we plot kW{x, k) at different locations, x, relative to 
the subshock, as indicated. The seed turbulence from 
the unshocked interstellar medium is at the level of the 
horizontal lines. 

Without cascading (solid curves), we see fairly narrow 
peaks centered at distinct k. This is in stark contrast 
to the case with cascading (dotted curves), and to the 
broad spectrum of k Wjx.k) seen previously with reso- 
nant interactions (e.g.. lVladimirov et al.ll2006D . Only the 
most energetic CRs reach a far upstream position near 
Xfob, and the small current these CRs produce gives a 
maximum 7nr at small k — kc/'^- At Xf^h = — lO^'r^o, 
where the light-weight curves are calculated in Fig. [3l 
7nr is a maximum at fc ~ Ix 10"'* r~Q with no cascading, 

and at fc ~ 5 X 10^"' r^Q* with cascading. The difference 
corresponds to the difference in escaping CR flux at a;feb- 
At a position within the precursor at x = — 1.2x lO^rgo 
(medium- weight curves), a second peak at larger k has 
developed when cascading is suppressed. At this posi- 
tion closer to the subshock, lower energy particles ap- 
pear in large numbers and the generation of waves at 
k ~ 10~^ Tgo shuts down because of the range in k given 
in Eq.[2l The increased number of lower energy particles, 
and the accordingly increased jd, now corresponds to a 
greater kc and shorter wavelength structures get ampli- 
fied at fc ~ 10^^ r'gQ . The position of the leftmost peak 

is nearly constant between Xi^h and — IxlO^rgo because 
id is nearly constant in this range, as indicated in Fig. [TJ 
At the downstream position (heaviest solid curve) , three 
distinct peaks are generated without cascading. 

The same turbulence generation model with cascading 
shifts energy from small to large k and produces smooth 
spectra. Far enough upstream in the precursor but not at 
Xfob (medium- weight dotted curve) , where only relatively 
long wavelengths arc amplified, the spectrum exhibits 
the "injection range" at fc < 10"'^ r~Q^, the broad "iner- 



tial range" 10 ^t^q 

spectrum W oc k~^/^, and the "dissipation range" at 
k > lOrgo- Closer to the subshock, the increased jd 
shifts turbulence amplification to greater k making the 
spectrum harder. The downstream spectrum becomes 
close to oc k~^ over a broad range in k. 

The peaks occur because of the coupling of parti- 
cle transport with turbulence amplification. The first 
(smallest k) peak forms far upstream, where only the 
highest energy particles are present, and their current jd 
is low. These particles diffuse in the X (x regime, scat- 
tered by the short-scale magnetic field fiuctuations that 
they themselves generate. As the plasma moves toward 
the subshock, advecting the turbulence with it, lower en- 
ergy particles appear. At some x, particles with energies 
low enough to resonate with the turbulence generated 
farther upstream (in the lowest k peak) become dom- 
inant. This strong resonant scattering leads to a high 
gradient of (seen at x ~ —10^ Vgo in the third panel of 
Fig- HI) J and the wavenumber kc/2, at which the ampli- 
fication rate 7nr has a maximum, increases rapidly. The 
increased value of fee /2 leads to the emergence of the sec- 
ond peak between 10~^ and 10~* r~Q , as seen in Fig. [3l 
Similarly, the third peak is generated at distances closer 
to the subshock than ~ — 1 x lO^'rgo and this is seen 
in the heavy-weight, downstream spectrum in Fig. [3] at 

The number of peaks depends on the dynamic range, 
i.e., on Xfeb- A smaller ccfeb can result in two peaks, while 
a larger a:feb, and therefore a larger Pmax, can yield four 
or more peaks in the downstream region. As mentioned 
above, Xfcb is a parameter in our model and is chosen 
here to give £;,„ax 10 - 100 TeV. 

The formation of the spectrum with discrete peaks oc- 
curs simultaneously with the stratification of the shock 
precursor into layers (see the plots of jd), in which vor- 
tices of different scales are formed. The peaks are a direct 
result of Bell's nonresonant instability, but they will not 
show up unless A and j^,- are calculated consistently, and 
the simulation has a large enough dynamic range in both 
k and p. 

As emphasized above, cascading in strong turbulence 
is uncertain and our two models bound the extremes of 
no cascading and Kolmogorov cascading to larger k. Our 
results are also limited by the fact that we haven't con- 
sidered the case of energy transfer from short to long 
scales. Nevertheless, the stark differences we see in both 
the shock structure and the turbulence spectrum suggest 
that observations may be able to constrain cascading 
models. The stratification is eliminated and the peaks 
are spread out if rapid Kolmogorov cascading is assumed. 
Downstream, the spectrum becomes close to, but slightly 
flatter than W oc k~^. Furthermore, the amplified effec- 
tive magnetic field, the shock compression ratio, and the 
maximum energy of accelerated particles are all signifi- 
cantly smaller with cascading. 

The large amount of energy observed in the first (small- 
est k) peak of the turbulence spectrum without cascad- 
ing defines a potentially observable spatial scale. For 

The tail at fc > 100 r~Q^ is produced by the large j^i at —1 Vgo < 
X < seen in Fig. [T] The short-scale turbulence produced by this 
current contains little power and does not impact our results for 
larger scales. 
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the parameters used here, -Dmax ~ 2-K/k « 2x lO'^r^o ~ 
10^^ cm, but this depends on Xfob, and we have found 
(work in preparation) that I?max oc \xich\ <x. pmax- 
Therefore, spatiall y resolved iiitensity variations trans- 
late to Umax (see iBvkov et alJl2008l i and may offer a 
new way to estimate Pmax from X-ray synchrotron ob- 
servations of SNRs. For our parameters, the turbu- 
lence on the scale Z^max has the effective magnetic field 
A_B « V AnkW « 400 /iG, and the corresponding Alfven 
speed VA ~ 400 km s~^, giving a characteristic time 
scale t = Dam^/vA ~ 0.5 yr. A larger Pmax implies a 
larger t. This rapidly varying field may influence both 
the spectral shape and the time variabil ity of X-ray syn- 
chrotron emission in S NRs (see, e.g., lUchivama et al.l 
[2007t lBvkov"eran[200l . 

On the other hand, if cascading to short scales is impor- 
tant, magnetic energy is efficiently transfered to the back- 



ground gas and the amplified magnetic field is consider- 
ably reduced. The heating of the shock precursor by the 
dissipated turbulence will also be significant. As shown 
in the bottom panel of Fig. [1] the precursor temperature 
can be orders of magnitude greater with cascading than 
without. The low-energy superthermal particles (up to 
p « Tripc) are more abundant in the presence of cascading 



and dissipation (to p panel of Figure 



, > P 

boosted injection (|Vladimirov et all 



H) du e to heating- 
2008f ). which may 



influence bremsstrahlung and spectral line emission in 
SNRs. 
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